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Abstract 

In this paper we introduce simplified, combinatorially exact formulas that arise in the vortex interaction 
model found in [29]. These combinatorial formulas allow for the efficient implementation and development 
of a new multi-moment vortex method (MMVM) using a Hermite expansion to simulate 2D vorticity. The 
method naturally allows the particles to deform and become highly anisotropic as they evolve without the 
added cost of computing the non-local Biot-Savart integral. We present three examples using MMVM. We 
first focus our attention on the implementation of a single particle, large number of Hermite moments case, 
in the context of quadrupole perturbations of the Lamb-Oscen vortex. At smaller perturbation values, we 
show the method captures the shear diffusion mechanism and the rapid relaxation (on Re 1 / 3 time scale) to 
an axisymmetric state. We then present two more examples of the full multi-moment vortex method and 
discuss the results in the context of classic vortex methods. We perform spatial convergence studies of the 
single-particle method and show that the method exhibits exponential convergence. Lastly, we numerically 
investigate the spatial accuracy improvement from the inclusion of higher Hermite moments in the full 
MMVM. 

Keywords: vortex methods, particle methods, vortex dynamics, Hermite functions, Navier-Stokes 
equations 



1. Introduction 

A new viscous vortex interaction model was introduced in [2"9l 157] whose main purpose is to study the 
dynamics of the viscous n- vortex problem. As a first application, the model was shown in [28] . to improve 
the approximation of the far field acoustic pressure field generated by a pair of rotating vortices, when 
compared to classical point vortex methods. One of the strengths of the model is the reduction of the 2D 
vorticity equation to a system of ODE's with purely quadratic nonlinearity which represent the evolution 
of the Hermite moments of each vortex. The main drawback of the model as presented in [22] is that the 
coefficients in the ODE's are computed in terms of derivatives and limits of an explicit kernel function which 
are computationally intractable for higher resolution (more Hermite moments) computations. 

We introduce in |Appcndix A| simplified, combinatorial exact formulas for the coefficients of these terms, 
greatly reducing the computational costs of generating and implementing the ODEs. These simplified for- 
mulas allow the viscous vortex interaction model of [29] to be implemented as a new multi-moment vortex 
method to simulate the two-dimensional vorticity equation. The novelty of this multi-moment vortex method 
(MMVM) is that the vortex particles are allowed to dynamically deform under the flow, without the usual 
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difficulties of computing the Biot-Savart kernel for anisotropic basis elements, thus reducing the convection 
error associated to classic vortex methods. 

It is important to note that the multi-moment vortex method we propose here shares many features with 
traditional vortex methods for incompressible [2l[7llllj and compressible [13 flows which are comprehensively 
reviewed in [T2] . One of the advantages that this method shares with traditional vortex methods is that it is 
gridless. In particular, we can avoid the common computational issues of using a large grid approximation for 
an unbounded domain. Using vortex methods, one also avoids the necessity of imposing either artificial, or 
periodic boundary conditions which may introduce unwanted or unphysical effects. MMVM also naturally 
incorporates the effects of viscosity in a way most closely related to the core spreading method, utilized 
in the work of Leonard, Rossi, Barba, Huang [331 [3J [T5]. Core spreading, using a purely Lagrangian 
vortex method, was originally shown to converge to the wrong solution in |15j . An adaptive method, 
proposed by Rossi [32| . has shown to fix this convergence problem and since then, innovative interpolation 
and initialization methods [5], have been developed to more accurately address this convergence problem. 

To incorporate anisotropic deformations, our method uses a Hermite moment expansion for each vortex 
element. For two-dimensional inviscid fluids, Melander et al. [33J HI] a ls° developed a model based on 
the evolution of moments of vortices. Instead of using Hermite moments, Melander et al. define a local 
coordinate system at each vortex and compute equations for the evolution of the centroid along with the 
time evolution of the local geometric moments. While in theory equations may be computed to any order 
for these expansions relatively low order approximations are used in |23[ 124] and convergence issues are not 
addressed. Using a Hermite moment expansion is natural, as they form a suitable basis and conditions 
for global in time convergence have been proven for MMVM [29 . In addition, it is assumed in (23J [H] 
that the maximum diameter of any localized vortex is much smaller than the minimum distance between 
any two vortices. This would specifically not allow the use of many overlapping vortices to simulate the 
fluid dynamics, i.e. a vortex method approach. Nonetheless, Melander et al. are able to show, using well 
separated vortices with moments up to second order, significant improvement over the point vortex method 
for several computed examples involving well-separated, co-rotating vortices. 

To our knowledge, the most closely related work to MMVM is Rossi's deformable vortex method using 
elliptically deforming Gaussian vortex particles. Rossi implemented the method first for the advection 
diffusion equation [33, 34] and also for the Navier-Stokes equations [35l[3T]. For both the advection diffusion 
equation and the Navier-Stokes equations, Rossi's deformable vortex method was able to achieve higher 
order convergence, as compared to non-deforming vortices, but only if one also modifies the Lagrangian 
velocity field. Unfortunately, the method was hampered by the difficulty of computing the Biot-Savart law 
for anisotropic (elliptical) basis functions to reconstruct the velocity field. While several ideas were proposed 
to alleviate this problem [351 131 j . we note that one of the most promising advantages of our approach is 
the avoidance of this difficult computation altogether, as the Biot-Savart integral is explicitly known to all 
orders in our Hermite expansion. 

For general flows the primary source of error from vortex methods arises from convection error which is 
caused by advecting vortices without deformation. Thus our method should greatly reduce this error since 
it allows for many different deformations, limited only by the order of the Hermite approximation. In fact, 
constructing basis functions that are allowed to deform have been conjectured as precisely what is need to 
improve the spatial accuracy [6] of vortex methods as seen already by Rossi's improvement from second 
to fourth order with elliptical deformed Gaussian basis functions [55] , Prior to Rossi's work elliptically 
deforming basis elements were used by Moeleker and Leonard [27 in a filtered advection-diffusion equation. 
In |27| the authors use a Hermite expansion but only to compute an approximation to the average velocity 
field, and not the vorticity field itself. This modification to the Lagrangian velocity field did not result in an 
increase the order of accuracy. 

In this paper we will first analyze the single-particle MMVM (n — 1) with a large number of Hermite 
moments, m. Unlike usual particle methods, our single vortex element must be large and model an entire 
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vortex structure. As a first application of the single-particle MMVM, we consider the problem of shear 
diffusion around a rotationally symmetric vortex. The shear diffusion manifests itself in this context through 
the mixing of vorticity and the axisymmetrization of the vorticity profile on a Re 1 / 3 time scale [8]. MMVM is 
a natural choice to study this problem for several reasons. As we will discuss in Section [2] our single-particle 
MMVM is equivalent to implementing a Hermite spectral method. Thus, while offering the usual good 
approximation properties of classical spectral methods it avoids the artificial periodic boundary conditions 
associated to classical spectral methods (Fourier, Chebyshev), see [5]. Aside from these general advantages, 
the single-particle MMVM offers some particular advantages for the study of two-dimensional fluid motion. 
Theoretical results have shown that the Hermite functions we use for our expansion give a natural basis 
with respect to which one can study the long-time asymptotics of solutions of the two-dimensional vorticity 
equation 14J. Furthermore, the vorticity distribution (in contrast to the fluid velocity) has the property 
that if it is initially localized it will remain so for later times. This means that if our Hermite expansion of 
the initial vorticity distribution converges, it will do so for all later times [29] , 

We next consider two examples using the full MMVM: A model for early time vortex merger behavior using 
n = 2 vortex elements, and second, we consider tripole relaxation using a course grid approximation n = 36. 
In both cases we observe increasingly accurate solutions as we increase the number of moments m. MMVM 
has several important parameters to choose for implementation. Just as in regular vortex methods, one 
must choose the number of particles n, and the overlap ratio which is computed by dividing the particle core 
size by the particle spacing. In addition, we must also specify the number of Hermite moments m for each 
particle to include in the simulation. As we will sec both parameters n and m can be increased to achieve 
greater accuracy and, at least in the case of a single particle (n = 1) the spatial convergence is faster than 
polynomial in m. 

We will work primarily with the 2D vorticity equation which can be derived from the Navier-Stokes equations 
for two dimensional incompressible fluids, 

| + (u.V)u = -V^+,Au, (1) 

V ■ u = 0, (2) 
where, u is the fluid velocity, p is the fluid density and v is the kinematic viscosity. 

If we take the curl of ([I]) and assume that the vorticity field u) = V x u is sufficiently localized then the 
equations for viscous vorticity on the entire plane are 

— — = vAui — u • Vw. 
dt 

uj = w(x,t), xel 2 , t > (3) 
w(x,Q) = wo(x) , 

where we can recover the velocity of the fluid via the Biot-Savart law 

u(x) = ±- [ (x ~ y ^ (y)rfy , (4) 
2irJ R 2 |x-y| 2 

and for a two-vector z = (z 1; z 2 ), z = (— z 2 , Z\). 

The rest of this paper is outlined as follow: In Section [2j we will briefly review the model from j^Hj which 
naturally defines the MMVM. In this section we begin with a review of the single-particle MMVM and then 
present the full MMVM. 

In Section[3j we present three different examples of our method. We first carry out the single-particle MMVM 
implementation in the context of perturbations of a monopole. We take our monopole to be the Lamb-Oseen 
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vortex and perturb the vortex by a quadrupole perturbation. We show, using an enstropy calculation, that 
the rapid shear diffusion mechanism causes axisymmetric relaxation to occur on the time scale of Re 1 / 3 , 
much faster than the purely diffusive time scale of Re. 

Our next two examples in Section [3] use the full MMVM. The first is an n — 2 model of vortex merger. By 
using two vortex elements to represent the merging vortices we show that increasing the number of Hermite 
moments improves the short time fluid dynamic behavior of each vortex. The third example is a coarse grid 
approximation (6 x 6) of a large quadrapole perturbation to the Lamb-Oseen vortex and we show that just 
including the second order Hermite moments significantly improves the early time simulation as compared 
to classic vortex methods. This example in particular illustrates the promise of implementing MMVM on a 
much larger scale. 

In Section|4j we perform two convergence studies on the single-particle MMVM and show that the convergence 
is exponential. We also quantify the improvement in spatial accuracy in the course grid calculation and end 
with a brief discussion of our computational implementation of MMVM. 

We conclude in Section [5] with a discussion of our results and future work in the further development of 
MMVM. Lastly, the simplified system of equations for MMVM is derived in |Appcndix A| 



2. Review of the viscous vortex model 



2.1. Single-Particle MMVM 



In this section, we review the viscous vortex model, presented in |29j . which expands the vorticity in terms 
of Hermite functions. We then formally state the single-particle MMVM, derived from the model in [29] . 
We begin with a brief review of the Hermite functions. 



Define 

^oo(x,t;A) = 4i e_|x|2/A2 ( 5 ) 

7TA Z 

where A 2 := Aq + Aid and Ao represents the initial core size of our localized vortex structure. Note that for 
any value of Ao, ^ defines an exact solution of the two dimensional vorticity equation known as the Lamb- 
Oseen vortex. As a consequence, we can choose any value of Ao in the definition of our Hermite spectral 
method; Ao is most often chosen to represent a typical length scale in the initial vorticity distribution. 

We now define the Hermite function of order (fc 1; /c 2 ) by 

4> klM (x, t; A) = D£i?£0oo(x, i; A) (6) 

and the corresponding moment expansion of a function by 

00 

w(x,t)= J2 M [fei> fe 2;*]^x,fc 2 ( x >*;A), (7) 

fei,fc 2 =0 

We define the Hermite polynomials via their generating function: 



H nun2 (z;\)= (D%D%e 



(8) 
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Note that the "standard" Hermite polynomials correspond to the case A = 1. Then using the standard 
orthogonality relationship for the Hermite polynomials, 

/ H til ,„ 2 (z;l)i/„ ll . m2 (z;l) e -^dz-7r2" 1 +" 2 (n 1 !)(n 2 !) ( 5„ l!ml( 5„ 2im2 , (9) 

JR 2 

we see that the coefficients in the expansion ([7| are defined by the projection operators: 

(•_l)fei + fe2^2(fei + fe 2 ) r 

M[k u k 2 ;t] = P kl ,Mt)] = 2 fci+fc2( fcl !)(fc 2 !) J R2 H *uk2(z;X)"(z,t)te ■ (10) 



If the function w(x, t) in Q is the vorticity field of some fluid, the linearity of the Biot-Savart law implies 
that we can expand the associated velocity field as: 

oo 

u(x,t)= M[k 1 ,k 2 ;t}V kuk2 {x,t;\) (11) 

where 

V klM (x, t; A) = d££^V 00 (x, t; A) (12) 

and Voo(x,i;A) is the velocity field associated with the Gaussian vorticity distribution </>oo- Explicitly we 
have: 

V M (x,,;A) = i(=5p)(l-e-M^) . (13 ) 

Remark 2.1. For simplicity, in this review we do not allow the centroid of vorticity to evolve in the single- 
particle MMVM. Of course, this is not the case in the full multi-moment vortex method. This restriction 
essentially adds the condition that the center of our expansion is chosen in such a way that both first moments 
are zero. 



2.1.1. Convergence 

One of the main theoretical results of [22] is the derivation of a criterion on the initial data which ensures 
convergence of the expansion ^ for all t > under the flow of 2D vorticity We state the proposition 
here: 

Theorem 2.2. Define the weighted enstrophy function 

£(*) = / oo 1 ( x ^;A)(^(x,i)) 2 dx . 

Jr 2 

If the initial vorticity distribution luq is such that £(0) < oo for some \q, and ujq is bounded (in the L°° 
norm) then the expansion Q is convergent for the initial value problem ^ for all times t > 0. 



The proof relies on a differential inequality argument and the fact that the Hermite functions are eigen- 
functions for a related self-adjoint linear operator, see for details of the proof. With this criterion for 
convergence, we now turn our attention to the governing equations for the moments which must be solved 
to implement the MMVM. 



2.1.2. Differential equations for the moments 



Assuming that the function w(z, t) is a solution of we note here that by using ansatz ([7 
clear that the single-particle MMVM is equivalent to a Hermite spectral method for solving 



for oj(z,t) it is 
31) , as the entire 
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effort is to compute the evolution of the coefficients of the Hermite functions. To begin we first introduce 
the following summation notation: 



J2 f(h,k 2 ):=J2 E f(h,k 2 ). 



ki ,k 2 



2 — fci+/C2— i 

/ci>0,fc 2 >0 



Thus we now define 



and 



l (x,f) = M[k ll k 2 ;t]4> klM {x,t;X), 



ki ,k 2 



u m (x,t) = J2 W[*i > *2;t]V fcllfca (x,t;A). 



(14) 



(15) 



(16) 



ki,k 2 



Here, uj" 1 and u m represent our Hermite approximation to solutions of Q. 

We can now derive equations for the evolution of the coefficients in the moment expansion by a standard 
Galerkin approximation. For nonlinear PDEs, this derivation can be challenging since passing the nonlin- 
earity through the projection operator ( 10 1 often cannot be explicitly computed. This is even more true 
in the case of vorticity where the nonlinearity includes a product of vorticity and velocity. However, the 
expressions for the coefficients in these expansions are surprisingly simple and explicit. Let P m [-] represent 
the projection on the subspace spanned by the Hermite functions of order m or less. (Without detailing the 
argument here, we consider this as a subspace of the weighted L 2 space in which the norm is defined by 
the weighted enstrophy function in Theorem 2.2 See [29] for a detailed discussion.) The standard Galerkin 
approximation to ^ is then given by: 



d t uj r ' 



E 



dM[ki,k 2 ;t] 
df 



tl! fc 2 (x,i; A) + ^2 M[k 1 ,k 2 ;t}d t (f>k 1 ,k 2 {x,t;X) 



k-\ ,k 



J2 M[k x ,k 2 ;t] {vAcf> klM (x,t;X)) 

m 

£Af[*i,4,;t]V/ lA (x,i;A) 1 V I E M[k u k 2 ;t](f> kuk2 (x,t;\) 



(17) 



_pr, 



i fcl ,k 2 



Our first simplification comes from noticing that the final term in the first line cancels the middle line of 
equation (17) (from the definition of Hermite functions) and hence if we apply the projection operators 
Pk lt k 2 , defined in (10 1, for k\ + k 2 < M we are left with the system of ordinary differential equations for the 
moments 



dM 



[h,k 2 ;t] 



-Pki,k 2 



^M[4,4;t]V^ 2 (x,*;A) 

/ m 

V M[mi,m 2 ;i\<p 



which, evaluated, takes the form, 



m m 



dM (_i)(ki+k 2 ) x2 (k 1+ k 2 ) 
^iM=- K - 2^(k^)(W) E E M[i u t 2 ,t]M[m M t] 

x / H klM (x)(D^D^V 0Q (x; A)) • V x (D^ 2 <Mx; A))dx. 

Js. 2 



(18) 



(19) 



G 



Explicitly computing equation |18|) is crucial to numerically implementing our single-particle MMVM (or 
equivalently our Hermite spectral method), since the projection operator is an integral over the entire plane. 
The difficulty in obtaining an analytic expression for the RHS of ( |18[ ) is that products of Vt u t 3 and <f>ki,k 2 
cannot be integrated in a straightforward manner. It was shown in |29] that this difficulty can be surmounted 
and equation ( 18 1 can be reduced to taking limits of derivatives of a relatively simple kernel function, K, 
irm: 

[ki,k 2 ,t] = 



thus taking the form 
AM 



where 



and 



dt 



(_l)(fcl+fe2)^2(fci+fc 2 ) 

2 fc i+ fe 2(fc 1 !)(fc 2 ! 



m rn 



J2 r[fa,k2jij2,m 1 ,m 2 ;\}M[£ 1 ,l2,t]M[m 1 ,m 2 ,t} 



<i,<! m 1 ,m 2 



T[ki, k2,h,£2,mi,m 2 ;X] = 

= D k T l D% D™* D{\ D e a lK(a 1 ,a 2 ,b 1 ,b 2 ,T 1 ,r 2 ;X)\ T =o.a=o.b=o 



K(ai,a 2 ,bi,b 2 ,ti,t 2 ;X) = -^ e -^ tiai+taaa) x 



1 

ttX 2 ' 

-a 2 fi + b 2 1\ + oi t 2 - bi t 2 
((&! + (h - fll ))2 + (b 2 + (t 2 - a 2 )f) 



1 — e 2a 



^2 ((fcl+(tl-ai)) 2 + (6 2 + (t2-a 2 )) 2 ) 



(20) 



(21) 
(22) 



(23) 



Implementing the single-particle MMVM is now reduced to solving equations (20 (-( 23 1 subject to appropriate 
initial conditions. Unfortunately, evaluating equation (21 1 directly is computationally intractable for large 
to. One of the main results of this paper is an explicit combinatorial formula for T[ki, k 2 ,£i,£ 2 , mi,m 2 ; A], 
derived in |Appendix A| Using the result from |Appcndix A| we can instead implement our single-particle 
MMVM by solving the system of quadratic ODEs given by equations ( A. 19 1-( A. 21 ). We will implement 
this single-particle MMVM in Section 3.1.1 to study the shear diffusion mechanism but next, we review the 
general MMVM model. 



2.2. The general MMVM 



To extend the single-particle MMVM we separate the solution of the vorticity equation into n vortex particles 
and derive separate evolution equations for each particle. Returning to the initial value problem ([3]), we begin 
by decomposing the vorticity distribution, writing (for t > 0), 

n 

w(x,t)=5^t^(x-x»'(t) > t), (24) 
j=i 

with velocity field 

n 

u(x,t) = ]T^(x-x^),t), (25) 

3=1 

where u J (y,£) is the velocity field associated to w J (y,£), centered at x J '(t), determined by the Biot-Savart 
Law. 

Remark 2.3. This decomposition is not unique. We can choose any number of pieces, n, into which we 
decompose the vorticity. However we will require two conditions on the choice of decomposition: the first is 
that the total vorticity of each vortex is non-zero, i.e. 

m? = I wg(x)dx^0, j = l,...,N. (26) 
Jr 2 
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The second condition requires that for t > 

/ (x-x J (t)V'(x-x J (f),i)dx = for all t>0,j = l,...N 
Jm. 2 



(27) 



As we will see, imposing condition \2D defines the motion of each vortex and is a departure from the 
traditional Lagrangian perspective of assigning the condition: 



dxj 
At 



u(x, t). 



(28) 



Physically, condition {2D imposes that the motion of the centers of the vortices are determined by momentum 
conservation, not by the passive advection of the local velocity field. 



There are two reasons why we favor ( 27 ) instead of the more traditional choice ( 28 ) . First, in \1J$ , it was 
proven that for the long-time asymptotic behavior of solutions of the two-dimensional vorticity equation, one 
has more rapid convergence to the Lamb-Oseen vortex solution if the initial data satisfies ( |27| . Secondly, 
support for this point of view is provided by the work of Lingevitch and Bernoff in \2(fj who argue that 
condition (27) is the "natural" definition for the center of a vortex immersed in a background flow. 



If we take the partial derivative of (24 1 and use the vorticity equation ^ , we find: 

n n 

d t uj(x,t) = ^<9 t ^(x-x J (<),t)-^x J (t)-Va; J (x-x J (t),i) 



i=i 



n / n 



J2 v&j (x - x? (*) , t) - J2 E u ' ( x - x " (*) > *) • Vuj3 ( x - xJ (*) > *) 

j=l j=l \Z=1 J 



(29) 



Given this equation, it is natural to define aj J as the solution of the equation: 



dt 



■{x-x 3 (t),t) = vAuj j {x-x j (t),t)- l^2v?{x-x?(t),t)\ ■Vu i (x.-x?(t),t) 



sl=l 



+*?(t)'Vu i (x-x>(t),t) , j = l,...,N. 



(30) 



It remains to determine x- 7 (t) . If we differentiate condition ( 27 ) we find that the equations of motions for 
the centers are defined as: 



dx» 
dt 



(t) = — V / «(z + x-''(t)-x £ (t),tV"(z,i))dz, n=l,2, 



(31) 



Equation ( 30 1 gives a set of n partial differential equations which govern the evolution and deformation 
of the vorticity of each localized vortex structure and equation (31 1 gives a set of 2n ODEs governing the 



motion of the centers of each localized vortex structure. Taken together, equations (30) and (31) arc the 



model proposed in [29 and solutions of the model solve the 2D vorticity equation exactly. 



2.2.1. The expansion for several vortex centers 



Let us now briefly comment on the extension of the Hermite moment expansion of the previous section to 
the case in which there are two or more centers of vorticity by combining this expansion with the multi- 
vortex representation of defined in the previous section. In this case, we consider the equations (30) for the 
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evolution of each vortex and then expand each of the functions ui J in Hermite moments as in the previous 
section. Thus, we define 

oo 

<j{m,t) = J2 M j [k 1 ,k 2 ;t} ( l> kuk2 (z,P,\) (32) 

ki,k 2 

for j = 1, 2, n. We make a similar expansion for the velocity field in terms of the functions V^,^, and 
insert the expansions into (30). Letting z = x — x J (i) and recalling that dt<f>k lt ki — v^-tyktMi we obtain: 



dM j [k 1 ,k 2 ;t] 



dt 

--Pfei.fc 



(33) 



n oo 



53 53 AP"[4^ 2 ;i]V £lA (z + s JjS i;A) 

(OO 
53 M J [mi,m 2 ; (z,t;A) 
mi ,"12 
/ oo 

(M;A) 



where Sjj/ = x J (i) — x J (i). 



Each of the above projections is fully simplified in Appendix ~K Thus, implementing MMVM is now a 



matter of simulating equations ( A.35 )-( A.42 ). The simplified combinatorial formulae derived in Appendix 



|A| allow for the efficient implementation of MMVM using basis functions constructed out of Hermite moments. 
Remarkably the convergence for the multi- vortex model is also established in [2"5] . 



3. Implementation and results 



In this section we implement the MMVM in three examples, the first being an implementation of the single- 
particle MMVM in the context of studying shear-diffusion. The second example is implementation of the full 
MMVM using n = 2 particles to model the early stages of vortex merger. The last example builds toward a 
larger scale simulation of the full MMVM by using a coarse grid (6 x 6) approximation to a tripole relaxation 
example. In the last two examples we specifically consider the effect of increasing the number of moments 
m used in the MMVM. 



3.1. The single-particle MMVM 
3.1.1. Shear diffusion 

It has been observed that small perturbations of an axisymmetric vortex monopole rapidly relax back to 
an axisymmetric state. It is has been argued [HI [H] that the shear diffusion mechanism is the cause of 
the rapid homogenization of small perturbations on the fast time scale Re 1 / 3 . Shear diffusion or the rapid 
relaxation back to an axisymmetric state was first associated to the passive scalar problem on the time scale 
Pe 1 / 3 , where Pe is the analogous Peclet number, by Rhines and Young in [3D]. Unlike the passive advection 
diffusion equation, the vorticity is coupled to the streamfunction, and yet Lundgren argued in [21] that the 
same shear diffusion mechanism also homogenizes vorticity on a time scale of Re 1 / 3 . Bernoff and Lingevitch 
further refined this assertion in [S], by showing that any non-axisymmetric, non-translational (i.e. zero and 
first Fourier modes) perturbation will decay on the Re 1 ! 3 time scale. Moreover, the higher the Fourier mode 
of the perturbation, the more rapidly the the perturbation decays. 
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This shear diffusion mechanism is based on a study of the linearized problem and thus is not guaranteed to 
hold for larger perturbation values. Here, we study small perturbations of the Lamb-Oseen vortex with purely 
Fourier mode m = 2 perturbation and for low to medium Reynolds number (Re — 500-4000). Using our 
single-particle MMVM we are able to capture a rapid relaxation to an axisymmetric state that is consistent 
with the shear-diffusion time scale of Re 1 / 3 . 

To begin, our initial perturbation of the Lamb-Oseen vortex (monopole) is of the form: 

o/(x,0) = A| x |2 exp( _^! )cos(20) . (34) 

Perturbations of this exact form (but large values of S) were used in the study of the formation of tripoles 
by Rossi et al. in [3)6] and by Barba and Leonard in [3J. It is not hard to see that if we also choose the 



total circulation to be M[0,0] = 1, then the initial conditions with the perturbation (34) can be re-written 
entirely in our Hermite basis as 

w(x, 0) = fo (x, 0) + 4<5(0 2O - <fo 2 ). (35) 

We point out here that we have implicitly chosen the intrinsic core size Ao = 2 for both the vortex and 
the expansion which we will now fix for the remainder of this section. In addition we define the Reynolds 
number to be Re = M[0, 0]/u = 1/v. From this we see that the initial data can be exactly represented 
in our expansion as: 

M[0,0](0) = 1, M[2,0](0) = AS, M[0,2](0) = -AS. (36) 



3.1.2. Results 



We now use an order 24 expansion to evolve initial conditions (36) under our single-particle MMVM to 



simulate the shear diffusion mechanism and axisymmetrization of the vorticity. Small perturbations of the 
quadrupole moment of the form (36) create elliptical deformations of the Lamb-Oseen vortex as seen in 
Figure [T] below: 



Initial Vorticity: S=.1 




n 



U 



0.07 



0.06 



0.05 



0.04 



0.03 



0.02 



0.01 



Figure 1: Initial elliptical perturbation of the Lamb-Oseen vortex of the form u)(x., 0) = 0oo(x, 0) + 4<5(</>2o(x, 0) — </>02(x, 0)), 
for S = .1. 
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To capture and quantify axisymmetrization, we adopt the conventions of [36] and define the nonaxisymmetric 
enstrophy of the vortex to be 



E 



(w(x)- < w(|x|) >) 2 dx 



(37) 



where 



< u x >= 



1 

2^ 



2tt 



u)(x)d6. 



This E represent the L 2 norm of the nonaxisymmetric portion of the solution and it is this quantity which we 



demonstrate decays on the Re 1 / 3 time scale. We solve the single-particle MMVM equations ( A.19M A. 21 ) 
found m|Appendix~A| using the initial conditions Q with S = .1 for Re = 500, 800, 1000, 1500, 2000, 3000, 



Air 
1 



and 4000. In Figure pi we plot E against unsealed time and observe that each solution indeed goes through 
a rapid axisymmetrization. 




Figure 2: Plot of nonaxisy metric enstrophy (E) vs unsealed time t, showing relaxation of purely azimuthal perturbation (second 
mode)using m = 24 moments. 



Notice in addition that in Figure [2j there are small oscillations after the initial relaxation in the cases of 
Re = 3000 and 4000. These oscillations capture the exchange of energy between the zeroth mode and the 
higher order modes that occurs before the diffusive time scale as discussed in [3] . To confirm that the shear 
diffusion acts on the Re 1 / 2, time scale, we rescale time. Figure [3] provides clear evidence that the early rapid 
decay collapses on the Re 1 / 3 time scale, confirming that for early times, this is indeed the appropriate time 
scale for shear diffusion. Further note that in Figure |4j when one tries to capture shear diffusion on the 
diffusive time scale Re, no such collapse occurs. In fact, the decay fans, out signaling a poor choice of time 
scale for the early decay. 

To observe the actual vorticity distribution, in Figure [5] we reconstruct the 2D vorticity at t = 200 for 
Re = 500, 1000, 2000, and 4000. Notice that diffusion dominates and smooths the vorticity for the cases 
of Re — 500 and 1000 as expected. Notice that we also capture short spiral arms, characteristic of rapid 
axisymmetrization, in the higher Reynolds number cases of Re = 2000 and 4000. 



3.1.3. A note on larger perturbations 

In the case where 5 is large, relaxation back to a monopole is unlikely to occur. Instead, it has been shown 
in [BJ |31 [55] that the vorticity distribution relaxes back to a rotating tripole (at least on an intermediate 
time scale). Although larger 5 perturbations create four distinct regions of vorticity and thus may not be 
well suited to a single-particle MMVM, we nonetheless obtain excellent agreement to the work in [6] in both 
frequency of rotation and magnitude for the low Reynolds number calculation Re = 500 and 6 = .25. 
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200 



Figure 3: Plot of nonaxisymetric enstrophy (E) vs t/Re 1 / 3 , showing relaxation of purely azimuthal perturbation (second 
mode)using m = 24 moments. 




t/Re 

Figure 4: Plot of nonaxisymetric enstrophy (E) vs t/Re, showing relaxation of purely azimuthal perturbation (second mode) 
using m = 24 moments. 
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Figure [6] can be directly compared to Figure 6.8 in [BJ. It is important to note that Barba's calculations were 
done using a very high resolution an adaptive Lagrangian scheme with thousands of vortex elements. The 
main differences in the results occur at the center of the tripole. We do not find a smooth elliptical center 
vortex and instead have a slightly "blobby" approximation. For higher Reynolds numbers, we do not find 
good agreement with the work of [BJ. This is likely due to our inability to resolve very small scale behavior 
without retaining many more moments, see Section [4] for further discussion. 



3.2. Full MMVM Examples 



As we saw in Section 3.1.1 a large single vortex element with many Hermite moments can closely capture 
the shear diffusion mechanism. One difficulty the method encounters is in resolving the finest scales, which 
includes filamentation. By dividing the initial vorticity distribution into many smaller vortex "particles" and 
keeping perhaps only a few moments for each particle, one could reasonably hope to better resolve the vortex 
dynamics similarly to the approach of classical vortex methods. In this section we present two examples of 
the full MMVM, the first being a simple n = 2 model for vortex merger which we now discuss. 



3.2.1. Modeling vortex merger 



In this calculation, we use the MMVM to model a classic problem of vortex merger which has been well 
studied [T51 [TJ [55J [SBJ [17]. In this example we use n = 2 vortex elements, initialized with non-overlapping 
Gaussian profiles centered at (—1,0) and (1,0). Our goal of this example is to observe the effect on the 
dynamics of merger produced by including the Hermite moments of each vortex. We specifically focus our 
attention on the effect that the number of Hermite moments has on the motion of the ccntroid of each vortex. 



Whether two axisymmetric, equal, and like-signed vortices will merge is determined largely by the ratio 
between the core size of the vortices, a, and the distance between them, b. For vortices with Lamb-Oseen or 
Gaussian profiles, the critical ratio which determines whether the co-rotating vortices destabilize and begin 
to merge is in the range of a/b = 0.22-0.24, see [55] HI HO] ■ If the ratio of the vortices are below this range, 





Figure 5: A plot of vorticity distribution at time t = 200. From top to bottom and left to right the cases are Re = 
500, 1000, 2000, 4000. 
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Figure 6: Perturbation vorticity normalized by oj max (0) with Re = 500 and 5 = .25. The frequency of roation and distribution 
of perturbation vorticity has very good agreement with Figure 6.8 in [6], The most significant difference is the more "blobby" 
and less elliptical center for later times. 
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then they are in the well separated regime and the vortices will co-rotate, diffuse, and their Gaussian profiles 
may begin to elliptically deform as they co-rotate. On the other hand, if Gaussian vortices are initialized 
above this critical ratio, the vortices rapidly distort, begin to merge, and one might observe ejection of spiral 
arms of vorticity. 

Here, we consider two different initial conditions for vortex merger: the first is well above the critical ratio 
at a/b — 0.375 and second is well below the critical ratio at a/b = 0.125. Each vortex is initialized with 
total circulation of M[i, 0, 0] = 1 and v = 0.001, and thus our Reynolds number is Re = M[i, 0, Q]/v — 1000. 
These two examples allow us to study the effect that increasing the number of moments for each vortex has 
on the motion of the centroid of vorticity. To compare both examples we non-dimensionalize time using the 
turnover time t* = 2 J*b 2 ^ anc ^ hnplement both examples with Hermite moments of order to — 0, 2, 3,4,5, 
and 6, from t* = to t* — .152, to focus on the early time behavior of merger. In the first case (a/b — 0.375), 
we expect the vortices to rapidly begin to merge. In the left plot in Figure [7J we plot the distance from 
the origin to the centroid of the vortices. In the context of the merger regime we expect to observe a rapid 
decline of this distance. 

1 .02 1 , 1 q 



0.98 
0.96 - 
co 0.94 - 
11 0.92 
0.9 - 
0.88 
0.86 - 

0.84' 1 1 L 

0.05 0.1 0.15 

t* 

Figure 7: In this figure we plot the distance from the origin to the centroid of the vortices. On left a/b = .375 and on the right 
a/b = .125. For both plots the m = plot is ( ), the m = 2 plot is (+), the m = 3 plot is (*), the m = 4 plot is ( ), the m = 5 
plot is (□), and the m = 6 plot is (o 

)• 

Let us first consider the to = trajectory in Figure [7| This is the case of a pure Gaussian basis element, 
and the constant trajectory shows that vortex continues to rotate at a strict distance of 1 from the origin 
and this is clearly inaccurate for this regime. As to increases, we see that the distance to the origin rapidly 
begins to decline and for early times, the trajectory of the centroid of vorticity's early descent to the origin 
seems to stabilize in m to a limiting trajectory. In particular we can see that in the time the vortices rotate 
by .152 radians the distance of the centroid to the origin has declined by 20%. We include a plot of the 
vorticity at t* = .127 for each to in Figure [5] 

There are several points worthy of comment in Figure [8] First, notice that the top left plot is the case 
of pure Gaussian vortices (to = 0). As expected, they simply turn and diffuse without any deformation, 
entirely unaffected by being inside the vortex merger regime. The next two plots on the top line, m = 2 
and 3, show large deformations and merger as the low order approximations attempt to exhibit qualitatively 
correct behavior. Notice that only elongated deformations occur for to = 2, since shearing and filamentation 
are not yet captured by only retaining up to to = 2 moments. For to = 3, shearing and merging have 
both begun but are only qualitatively correct. By to = 4, the vorticity distribution is beginning to settle 
into the approximately correct vorticity distribution. Focusing on the to = 6 vorticity plot, we see all the 
characteristics expected in the early onset of merger: shearing and the early development of spiral arms of 
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vorticity, deformation of the initial circular core near the centroid of vorticity and mutual, rapid attraction 
of each centroid of vorticity toward the other. This example, while only a model for vortex merger, clearly 
demonstrates that it is precisely the inclusion of higher Hermite moments that allow us to more accurately 
capture the correct fluid dynamics of vortex merger. 



In our second example, we consider the case where the ratio of the co-rotating vortices is a/b = .125, well 
outside the critical ratio. Here, we expect that the centroid of vorticity should maintain the distance from 
the origin or perhaps only slightly decline, and this is precisely what we observe in the right plot in Figure 
[7j In addition we expect the well separated vortices to maintain there circular shape and perhaps begin to 
a exhibit slightly elliptical deformations. We plot the vorticity distribution at t* = .127 in Figure [9] First 
notice that in each plot of Figure [9] the rotation angle is the same. There is also very little deformation of 
the vortices with perhaps the slightest elliptical deformation being exhibited in the distributions with higher 
Hermite moments as expected. 



These two examples demonstrate that simply through the inclusion of higher order Hermite moments, our 
model can capture significantly more accurate physics. In the m = case of diffusing Gaussians, the 
governing equations for short times are essentially equivalent to the point vortex theory which states that 
two point vortices of equal circulation will rotate on a circle. Thus, it is precisely the inclusion of the higher 
order Hermite moments that allows two initially Gaussian vortex elements to exhibit convective merger. 
Moreover, both examples demonstrate that the Hermite moments correctly select the regime of the vortices 
(co-rotation or merger). Of course, vortex merger is an extremely complicated behavior and to continue 
to exhibit accurate solutions for longer times, more Hermite moments are ultimately needed to capture the 
dynamics. 



M=0 M=2 M=3 




Figure 8: The plots of the vorticity in the merger regime at t* = .127. The plots (from left to right and top to bottom) are 
m = 0,2,3,4,5,6. 
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3.2.2. Tripole Relaxation 



In our final example, we present a coarse grid direct numerical simulation of large quadrupolc perturbations 

As noted in Section 



of the Lamb-Oseen vortex of the form (36) found in Section 3.1.1 



3.1.1 larger values 



of S cause the vorticity distribution to relax to a rotating tripole, see |6j [3l [36] for careful studies of this 
phenomena. We have already shown in Figure [6] that a single vortex element with to — 22 moments 
successfully captures this relaxation for lower Reynolds number (Re — 500). We now revisit this problem, 
but at higher Reynolds number(i?e = 1000), and use a coarse grid (6 x 6) approximation of the initial 
vorticity distribution. For this example, we use an initial Lamb-Oseen vortex with core size Ao = 1 and our 
initial vorticity distribution once again takes the form: 



W (x) = w(x, 0) = 000 (X, 0) + 45(020 - 002)- 



(38) 



For this experiment, we select S — .25, a large quadrapole perturbation of the Lamb-Oseen vortex (note that 
this perturbation is even larger than in Section 3.1.3 since Ao is smaller). We overlay our 6x6 grid of vortex 
centers on [—1,1] x [—1,1] and approximate the initial condition (38) by initializing the circulation of each 
initial Gaussian vortex element using the standard approach 



M J [0,0;0] = wo(x i )Ax j 



(39) 



where Axj is the area of each node. We first compute a pure Gaussian simulation (to = 0), which is 
equivalent to a standard core-spreading vortex method, and then for our second simulation we re-run the 
identical initial conditions but include quadrupole moments (to = 2) for each vortex element. Since we are 
using core-spreading and a coarse grid, we focus our study on the early relaxation dynamics. A plot of the 



coarse grid approximation of (38) is presented in Figure 10 



In Figure [TT] we juxtapose the implementations of the full MMVM at t = 15 for to = and to = 2 using 
the same initial conditions |39[ The left plot is the pure Gaussian to = computation. In the center plot is 
the result of including the quadrupole Hermite moments to = 2 into the calculation. The right plot is high 
order (to = 24) single-particle MMVM simulation of the same initial conditions. 
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The results shown in Figure 11 arc quite striking. The right plot is taken to be a high accuracy solution 
of the identical initial conditions at t — 15 and captures the early time tripole relaxation observed in the 
work of Barba et al. O 0] and Rossi et al. [36] . In this plot we observe shearing and relaxation of the small 
satellite vortices. In addition, at each contour level we see characteristic spiral arm wind up occurring as 
well. 



With such a coarse grid, we expect to miss much of these important characteristics of the early tripole 
relaxation. Indeed, in the case of classical vortex methods with Gaussian basis functions, i.e. MMVM with 
m = 0, the approximation is poor. There is barely any spiral arm wind up occurring in the distribution of 
vorticity, and very little shearing and relaxation is observed. The simulation with m = looks more like a 
rotation and diffusion of the initial condition than the characteristic nonlinear rapid relaxation to a tripole. 

In contrast, for m — 2, the MMVM calculation shown in the middle plot of Figure [Tl] shows a dramatically 
more accurate approximation to the high accuracy solution on the right as compared to the m = example. 
In this case, we observe spiral wind up at each contour level of vorticity and the small satellites of negative 
vorticity have started to shear the outer arms of positive vorticity which is also characteristic of the early 
relaxation dynamics. Since these two experiments are initialized with identical initial vorticity, we can 
attribute the improvement of accuracy in the early relaxation dynamics of a tripole solely to the inclusion 
of the higher moments. We will quantify precisely this observed improvement of accuracy in Section [4j 
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Figure 10: The initial vorticity distribution constructed from equation |39| 
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Figure 11: The plot at t = 15 on the left is the solution using MMVM with m = and the plot on the right is the solution 
using MMVM with m = 2. 
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4. Convergence studies and performance 



In this section we will present two convergence studies for the single-particle MMVM: The first test case 
demonstrates the ability of our method to accurately approximate the Lamb-Oseen Vortex, an exact solution 
to ([3]) . The second example is a tripole relaxation example where no exact solution is known. In each example 
exponential convergence is demonstrated numerically. We will next quantify the improvement of accuracy 
in the coarse grid approximation example presented in Section |3.2.2| Finally we will briefly discuss the 
computational implementation used to compute these examples. 



4.1. The single-particle MMVM 



4.1.1. First test case: the Lamb-Oseen vortex 



A standard test example for 2D vortex methods is the Lamb-Oseen monopole which has initial conditions: 



w(r,0) 



-To - fl 
e 



4nX 2 



(40) 



where r = (x 2 + y 2 ) 2 . Initial conditions (40) are particularly useful because it results in the only known 
nontrivial, analytic solution to system ^ known as the Lamb-Oseen vortex. This solution takes the form 



v ' 4tt(\ 2 + Avt) 

and this example is often a first benchmark to compare the spatial accuracy of numerical solutions, [4j 135]. 



Provided that we select our core size Ao in our expansion to be equal to the initial core size of the initial 
data (40) then, by construction, the Lamb-Oseen monopole represents the zeroth order moment in our 
expansion M[0, 0], which is constant because the total vorticity is a conserved quantity [52]. Thus, no 
nontrivial dynamics will occur if we were to start a single-particle MMVM simulation using initial conditions 
M [0, 0] (0) = 1 and M[i,j] (0) = for all i, j; in this case, our method predicts that for all time t, M[0, 0] (t) = 1 
and M[i,j](t) = for all 7^ (0,0), in agreement with the exact solution. Hence, to perform a more 

meaningful convergence study for MMVM, we must approximate the Lamb-Oseen vortex in a way which tests 
the dynamics of the higher moments as well. We therefore choose to approximate a Lamb-Oseen vortex with 
core size Ao = 2.1, which remains an exact solution, with a Hermite expansion using a Hermite expansion 
with core size Ao = 2, thus allowing the higher moments to accommodate the error between the smaller core 
size Gaussian and the true solution. For our convergence study we select Re = 1000. 



We calculate the error between the true solution for the Lamb-Oseen monopole and the m moment approx- 
imation obtained with the single-particle MMVM for a range of values of m. In Figure [12J we show the 
L°° norm of the relative error plotted against m. The linear nature of these plots confirms that, as one 
expects from a spectral method, the error between the true solution and the simulated solution converges 
exponentially as m — > 00. The overlap of the four lines corresponding to times 2, 4, 8, 16 and 32 indicate 
that the error varies extremely slowly in time. This is likely due to the axisymmetry of the initial data and 
the fact that diffusion is the sole component of the dynamics. 



4-1-2. Test Case 2: A Larger Quadrupole Perturbation 

The previous test case numerically demonstrated that the convergence rate is exponential in the number of 
moments, in the case of axisymmetric initial data. We now test our method for non-axisymmetric initial 
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conditions which result in nontrivial dynamics for the evolution of the coefficients in our expansion. To this 
end, we select as a second test case, large quadrupole perturbations of the Lamb-Oseen vortex with initial 
conditions of the form (36) which were studied in Section 3.1.1 We select large enough perturbations such 



that the solutions do not relax to an axisymmetric state and instead evolves toward a metastable tripole 
solution, see [31136). which has no known analytic solution. Therefore, in order to measure the convergence of 
the single particle MMVM, we take a self consistent approach by selecting a high-order simulation (m = 24 
moments) in lieu of an exact solution and compare our lower-order simulations, in increasing order, against 
this solution. 



Log 1Q of sup norm of relative error plotted against m for early times 




Figure 12: The log 10 of the L°° norm of the relative error between the Lamb-Oseen monopole with core size 2.1 and our test 
approximation with An = 2 plotted against m, the number of moments retained in the simulation. The different lines in the 
figure correspond to times 2, 4, 8, 16 and 32. 




Figure 13: Here, log in of the L°° norm of the relative error of the solution is plotted against the number of moments m for 
Re = 1000 and <5 = 0.25. The different lines correspond to times 2, 4, 8, 16 and 32 (see legend). 
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In Figure 13 we examine the log 10 of the relative difference between the simulated solution with 24 moments 
and the simulations with fewer moments for 6 — 0.25 and Re = 1000. The different lines in the plot 
correspond to times t = 2, 4, 8, 16 and 32. The linear nature of the plot once again indicates that the 
L°° norm of the relative error decreases exponentially in to. The rate of exponential convergence is higher 
for earlier times. In addition, the error is very small at early times but increases over time. As expected, 
for any fixed number of moments, the dynamics will eventually diverge from the "true solution," however, 
it is clear from this figure that the rate at which this divergence occurs depends heavily on the number of 
moments we retain in the simulation. Thus, in practice, one can tune the number of moments kept in order 
to approximate the solution at the desired level of accuracy. 



In addition to examining the size of the error as to increases, it is interesting to consider the spatial distri- 
bution of the error. In Figure [hi] we plot in the icy-plane the relative difference, 

u n {x,t) - 0J2i{x,t) 
1 1^24 1 |oo 

between the solution for 24 moments and the solutions for a smaller number of moments at time t = 32. 
These numerical results suggest that in addition to increasing overall accuracy, the length scale of the error 
is declining as more moments are retained, indicating that resolution of the finer length scales improves as 
to increases. 



4-2. Error analysis of the coarse grid simulation 



In section 3.2.2 two coarse grid approximations of early time tripole relaxation were computed using MMVM, 
the first a was classical vortex method (using Gaussian basis functions without higher Hermite moment 
corrections), i.e. to = 0, and the second was a to = 2 MMVM simulation. When compared to both high 
precision tripole calculation, [BJ [3] and a high order (to = 24) single-particle MMVM calculation, the to = 2 
solution captured much more of the physical features associated to early time tripole relaxation as shown 
in figure 11 To quantify this improvement in error of to = 2, as compared to m = 0, we computed the 
relative L 2 , and L°° error, which we denote as ||e m ||2 and ||e m ||oo respectively for m = 0,2. We use our 
single-particle MMVM to = 24 calculation as the benchmark solution and present the result in Table [T] 
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0.1742 


0.0821 


0.1532 


0.0956 



Table 1: The relative error 6?-, 



and 



for m = 0, 2 at times t = 1, 2, 4, 8, and 16. 



It is immediately clear from the results in table [T] that the to = 2 MMVM simulation is a significant 
improvement in accuracy as compared to to = 0. Specifically we can see that for early times the to = 2 
simulation is approximately 4 times more accurate than to = in both L 2 and L°° metrics. This large 
improvement in accuracy eventually declines to around twice as accurate by t = 16. This decline in accuracy 
is due to the increased excitation of the second order moments which become too large and is analogous to 
the reduction in accuracy of Rossi's elliptical vortices when the aspect ratio becomes to large, [35]. Thus, 
like Rossi's elliptical vortices, to maintain higher order accuracy over longer time scales one will need to stop 
the calculation and re-interpolate or re-initialize the vorticity, likely using the high order accuracy methods 
recently developed in [5]. 
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Figure 14: The relative error at time 32 between the solutions for (from left to right) m = 4, m = 6, m = 8, m = 10, m = 12, 
m = 14, m = 16, m = 18, m = 20 and m = 22 and the solution for m = 24, with & = 0.25 and Re = 1000. 
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4-3. A note on implementation and performance 



In our implementation of MMVM, we symbolically compute the equations ( A. 19 1-( A. 21 ), in the case of the 
single-particle MMVM, or equations ( |A.35[ )-( [A42| ), in the case of the full MMVM, using Maple 10.5. The 
equations are then exported from Maple to Fortran code. We proceed to use this Fortran code to solve 
the equations in Matlab using ode45 with 'RelTol' and 'AbsTol' both set to 10 -8 . There is a one-time 
preprocessing cost associated with Maple's generation of the Fortran code and the subsequent writing to .m 
files. Once the .m files are created, the ode45 solver is relatively fast. For example, a typical run with initial 
conditions of the form (36) and Re = 1000 with m — 2 moments took less than one second, for m = 18 
moments took 3.35 minutes, and for m = 22 moments took 10.65 minutes to solve from t = to t = 2000 
on a 2.4 GHz Intel Core 2 Duo processor with 2 GB of RAM. Typical run times on the 6x6 course grid 
calculation with m = 2 took about 6 mins to run from t = to t = 50. These run times reflect the fact that 
there where several limitations of this simplistic implementation. Much improvement may be achieved as 
no effort was made to optimizing the memory usage and the efficiency of storing and calling the system of 
equations (which can become quite large) . New, efficient and parallelizable code is currently in development 
and is a necessary step for future large scale simulations using MMVM. 



5. Conclusion 



In this paper, we have presented simplified equations for a new multi-moment vortex method for computing 
solutions to the 2D vorticity equation. The novel feature of MMVM as compared to classical vortex methods 
is that the higher Hermite moments allow the vortex particles to convect and deform without any of the usual 
computational difficulties associated to calculating the Biot-Savart kernel for anisotropic vortex elements. 

We first considered a single-particle MMVM simulation and showed that the method captures the physical 
mechanism of shear diffusion (and the associated time scale Re 1 / 3 ) for low to medium Reynolds numbers. 
We then presented two examples using the full MMVM which show how the inclusion of higher moments 
improves the calculation over non-deforming basis elements. The first is a simple model for vortex merger 
and we show that early convective merger behavior is only captured with the inclusion of higher Hermite 
moments. In the second example we perform a coarse grid calculation of tripole relaxation and show that 
by including just m = 2 Hermite moments we show upwards of a four fold improvement in the reduction of 
error as compared to tradition vortex methods (m = 0). 

One interesting feature about MMVM is that there are two parameters, the number of particles, n, and the 
order of the Hermite expansion, m, that can be used to tune the spatial accuracy of the method, whereas 
classical vortex methods has just the number of particles n. It was shown in Section [4] that for a fixed n = 1 
the spatial convergence is exponential in m. One important avenue for future work is to better perform 
convergence studies and analysis of the spatial accuracy as a function of both n and m. For a fixed n, 
increasing m also introduces additional computational costs and in any practical implementation, one will 
need to find an optimal balance between computational efficiency and spatial accuracy by varying both n 
and m. 

The examples in this paper and the demonstrated improvement in the reduction of spatial error, establishes 
the promise of MMVM for larger scale calculations and thus the need for future work in this direction. In 
addition to an efficient implementation of MMVM, we will also need to incorporate spatial adaptation of 
the vortex particles, such as the highly accurate methods recently developed in [5], in order to push MMVM 
simulations to later times while maintaining accuracy. 
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Appendix A. Combinatorial formulae 



Appendix A.l. The Single Vortex Case 



In this appendix we calculate explicitly the integral term in equation (19 1 found in Section [2j which has the 
form 

f H klM (x)(D^D™ 2 V oa (x; A)) • Vx^^-Mx; A))dx . (A.l) 

Jr 2 



It is the goal of this appendix to further improve upon the simplification ( 23 ) by deriving exact combi- 



natorial expressions formulae for (A.l) since such combinatorial formulae will allow for efficient numerical 



implementation of the Hcrmite spectral method. 



We proceed by apply integrating by parts to ( A.l ) to get: 

/ iJ fel)fc2 (x)(£»^£)™/Voo(x; A)) ■ V^^^x; A))dx = (A.2) 
- I ((D^D™Voo(x;X)) ■ V x H kuk2 (x))Di\Dil4> 00 (x;X)dx = - (T l + T H ) (A.3) 

JR 2 

Remark Appendix A.l. The term in the integration by parts formula in which the derivatives fall on 
Vqo do not appear because the fluid velocity is divergence free. 



where 



T 1 - / (^ 1 1 ^V o 1 o(x;A))^ lJ ff fcl , fc2 (x)^ iJ D^0oo(x;A)dx (A.4) 

JR 2 

r" = / (^^V o 2 o (x;A))^ 2J ff fcl ^ 2 (x)^ i ^ 2 O o(x;A)dx (A.5) 
Jr. 2 

We first focus on equation T 1 . If we again integrate by parts l\ times with respect to X\ and 1% times with 
respect to x 2 we arrive at: 

rI = EE(t)(t) M) ' 1+£2 / ^ooD^Di 2 H kuk2 D^^D^ 2 -^doc (A.6) 
i=0 j=0 ^ / \j / Jv 2 

mirj(^i,fci-l)min(ij,fc2) , * , ni+li i \ / ni; t 



x 



b as H kl _ i _ 1M _ j D^- i D^+ t '-iV^dK 
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Where in the second equality we used the fact that 
Now if we take a closer look at the last integral we can compute: 



2n 

2m 
A2" 



/ j. tj rymi+i-L—i jymv+l*—j\r\ J„ 

/ Poo-Hki-i-lM-j-Vx! u x 2 v oa ayi 

JR 2 



JR 2 



jymi+ki-i-l+£i-i ]jm2+k2-j+i2— ^Y^dx 



where the first equality comes from the fact that 

H n , m = (-1)"+™^ D^D™ fa, (A.7) 
and the second inequality is a result of applying integration by parts. Thus in order to better understand 



expression ( A.l I we must simplify integrals of the form: 

/ ^)D^D a b 2 V^{x + b)dx| b=0 

JR 2 



(A. 



To do so we consider both components of the velocity field. By re-writing the velocity Vqo in terms of the 
vorticity 4>oq we have: 

V oo (x + b;A) = -V^(A h )- 1 0oo(x + b) , (A.9) 



where Vtf = (d X2 f,—d Xl f). Thus we can study equation (A.8) as the first component of the equation 
(lATOl below: 



KK [ V oo (x + b;A)0 oo (x;A)dx 

" JR 2 

= VJ(A 6 )- X / oo (x + b;A) ( / )oo (x;A)dx. 

JR 2 



(A.10) 



As seen in [29] we can reduce (A. 10 1 to the form 



f <p m D^D^V o(x + bA) = D%D%VooQ>; V2X)\ h=Q 

JR 2 



(A.11) 



Remark Appendix A. 2. It is worth noting here that the differential equations governing the moments of 
the Hermite representation reduce to evaluating derivatives of the velocity field associated to the Lamb-Oseen 
vortex at zero. As we will see in the next section similar importance on the velocity field of the Lamb-Oseen 
vortex persist in the multi-vortex expansion. 
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We consider the first component of (A. 11) which is relevant for equation (A. 8). Expanding this component 
we get 

2n b* ^ J 2n ;M2A 2 J (n+1)! 1 1 



n=0 



(A.12) 



Thus we can evaluate: 

1 (-b 2 ) f, „_M 
2tt |b| 



K^fe^Tuir 1 - e"^ | b=0 (A.13) 



_1 ~ / i (-1)" ^fn\ (2k)l 2k _ ai (2n-2fc + l)! , 2 „_ 2fe+1 _ Q2| 

^M^J (^l)!^ W (2fe- ai )! (6l) (2n-2fc + l- a2 )! (62) |b =°- 



n=0 v ' v ' fe=0 

We can then see that the relevant values of k* and n* are: 

k* = si 

2 

„ a 2 + an — 1 

n = . 

2 

It is also clear that a\ must be even and a 2 must be odd in order to have a non-zero contribution. This 
yields the final formula for the first component: 

U b 1 U b 2 7r-^rrr 1-e » | b =o - 



2tt |b| 

-1 / 1 \ n * +1 (-1)"* /n*\ (2ft*)! (2n*-2ft* + l) 



(A.14) 



2tt V2A 2 J (n* + 1)! / (2ft* - a x )! (2n* — 2ft* + 1 — a 2 )! 

1(01)1(02)! (A.15) 



-1 / 1 (-!)"' (n* 



or. 



where 



2tt V2A 2 J (n* + 1)! \k 



Q 2 + ai +l °2+°l-! Q2+Q1 _i 

Hi(ai,a 2 )=<| iKsW 2 ( °^ + 2 ° 1 + 1 )! ^ 4 K" 1 ^" 2 ) 1 if even and a 2 odd 

otherwise. 



If we compute the second component of (A. 10) in a similar fashion we see that 



D^D^^-^(l-e-^) b -„ = « 2 (n 1 .n>) (A. 1.7) 

where 



»2 + °l+l °2+gl-l f , a+ai -l 



= J ^(d?) 2 ( f 'Ui+i v ( -lLi )("i) ! ( Q! 2)! if a x odd and a 2 even 



otherwise. 
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We can plug these "simple" formulas into r 1 = T 1 [k\ , k 2 , i\ , i 2 , m\ , m 2 ; A] and T n = T 11 [ki , k 2 , l\ , l 2 , mi , m 2 ; A] 
to get our finalized equation for the coefficients of the moments, though we must first note that if we follow 
the derivation of I we can find that T n takes the form: 

r n [h,k 2 ,£ 1 ,£ 2 ,m 1 ,m 2 ;X] = £ £ ( h ) M (-l)^ 2 f ^^^H^D^+^D^^V^^x 

min(i lt kx) min(e 2 M 2 -l) /£ \ /£ \ / 2*fc ! \ / 

S 5 l^(j 2 j (_1) ' 1+ ' 2 (^w(fc^j ( 



2J' +1 A; 2 ! 



x / ^Hk^M-^D^^DZ^^V^ (A.18) 

JR 2 



Thus we can use the equations above to write for both T 1 and T 11 : 



T^M^mum^X] = E E Him H)<' ' — W — ^ 



i=0 j=0 



min(£i,ki — l) min(l 2 ,k 2 ) 



53 E ^ ^^(fci-i-i)!; u 2 M(jt 2 -j)! 



x 



i=0 j=0 

Hi (mi + fci - i - 1 + £1 - i, m 2 + k 2 - j + £ 2 - j) (A.19) 



min(li,ti)min(f 2 ,ti-l) . . . . . . . , . 



X 



i=0 i=0 



e e (XV»W 



2^1! \ / 2J +1 A;2! 



iJViJ \x 2 ( i ){ki-iy.J V^ 20+1) (fc2- i-i)! 



X 



i=0 j=0 

W 2 (mi + fci - i + 4 - i,m 2 + A; 2 - j - 1 + f 2 - j) (A.20) 

Thus we have simplified the differential equations for M[k\, k 2 ](t) to: 

— [ki,k 2 ,t] = 2 i 1+k2( , ,w, p. E E ^' m 2 ,t]r[fc l ,fc 2 ,4^2,m ll m 2 ; A], (A.21) 

where T[ki, k 2 , li,l 2 , mi, m 2 ; A] = r I [fei,fc 2 ,4,^ 2 ,mi,m 2 ;A] + r II [fci,fc 2 ,£i>^ 2 ,mi,m 2 ;A]. 

Appendix A. 2. The Multi-Vortex Expansion 

We continue our work by analyzing the multi- vortex version of the model in a similar way. First we re- write 



equation ( 30 ) which are the governing equations for the multi- vortex model: 



Ft i I m \ 

' W -(x-x^),t) = uAu} j (x-x^(t),t)- ^V(x-x*(i),i) -Vw^x-x^t),*) 



0t 



+x i (t)-Vw 5 '(x-x''(f),t) , j = l,...,M. (A.22) 
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Each of the functions w J is then expanded in Hermite moments as done in the previous section. Recall also 
that we must couple equations (A.22) with the equations of motion for the centers defined above in equation 
El) by: 



where Sjj> = — x J (t) + x J (i) and z = x — x J (i). 
To start we define 

oo 

u j {z,t) = M J '[fci,* 2 ;t]^ 1 ,fc a (a,t;A) 



(A.23) 



(A.24) 



for j = We make a similar expansion for the velocity field in terms of the functions V^^ 2 , and insert 

the expansions into (A.22). Continuing to follow the work in section 4 of [29] we find that the differential 



equation associated for each coefficient is: 

dMi[k 1: k 2 ;t} _ 
dt 

-PkiM 



' V ( X! M j [mi,m 2 ;t]cj) mitm2 (z,t;X) 

\mi ,m 2 
/ oo 

x J (t)-Vl ^2 [mi, ma; (z,t;A) 



(A.25) 



Equation (A.25) will be analyzed in several pieces. First note that when f = j then Sjji — and this 



case is reduced precisely to the work done in Section ??. Thus there are only two new pieces that must be 
computed, the vortex interaction terms which include in the first term the case j ^ j' and the entire second 



term in (A.25 ) 



For the case where j ^ j' , a similar derivation found in the previous section will retain the same conclusions 



up until one reaches the step at equation (A.ll). We set b = Sj t j> instead of b = here. This is as far as 
one can simplify the equations for the case j ^ j' . Thus we are left with the last term: 



x J 'W-V( M j [mi ) m 2 ;t](l> mi , m2 {z,,t;\) 

\mi,m 2 / 

Mi[m u m 2 ;t}[ H klM {z)V(t) ■ V z ( J D™ 1 D™^ o(z))dz 



(_l)(fcl+fc2) ^2(fei+fc 2 ) 



= III + IV 



where 



III = a^(t) / i/ fel! / £2 (z)Ar ll +i,m 2 0oo(z)dz 

/ -^fel ,k 2 ( z )-Dmi+l,m 2 0OO 

(z)dz. 



IV 



(A.26) 



(A.27) 
(A.28) 
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For now we consider just III. We see that by applying integration by parts we have: 

III = xi(t)(-l) m i +m * +1 [ J D mi + 1 '" i2 i7 fcl , fc2 (z)0 oo (z)dz 

Jr 2 

( lV n 1+ma + l ( 2 "' 1+lfc l ! A ( 2 "' 2fc 2 ! 

' I A 2(m 1+ i)( fcl _ TOl _ !)i J l v A 2 ("^)(fc 2 - ma)! 



x / H kl (z)^oo(z)dz. 

Jk 2 



But notice that 



/ 



, , , , i 1 ai, a2 = 
ff Ql , Q2 (z)<Mz)dz^ Q otherwige 



Similar work can be done for IV and thus getting for both: 



m= l - 1 (SS^)(^r)*i(*) *i = mi + l, fe = m 2 

otherwise 



(A.29) 
(A.30) 

(A.31) 



IV = 



■1) 



m 1 +m 2 + 



1 / 2 m ik 1 \ \ / 




2 ^(ll+ 1 ) ) fc l = TO 1> fc 2 = ™2 + 1 

otherwise 



Combining the equation (A.31 1 and (A. 32 1 with (A. 26) we get the remarkably simple equation 



(A.32) 



Pt 



fei ,fe 2 



x'W-Vl M 3 [mi,m 2 \t]<j> (z,i;A) 



i{(i)M J [fc! - l,fc 2j t] +4(t)M J '[/ci,fc 2 - l,i]. 



(A.33) 



Lastly we must analyze equation (A.23). Again we sub in our expansion (A.24) into (31) and find compo- 
nentwise: 



dx^ 
~dT 



(t) 



Mi [0 



1 m f 



z + Sjj/ , t)u^ (z, t) J dz 



m oo oo 



L ' ' ' j'=i e u e 2 mi,n2 

x / V£ 1 .£ 2 (z + s JJ /,i;A)0 mi . m2 (z,t;A)dz , 
Jr 2 



(A.34) 



where we write out 



f V, 1 ., 2 (z + b ! i;A)0 mi . m2 (z,i;A)dz = f rftofc V 00 (z + b)^^0 O o(z)dz 

JR 2 JR 2 

(-l) mi+m ' 2 / < +mi ^ +m2 V oo (z + b)0 oo (z)dz 

JR 2 



and now we are back precisely to equation ( A. 11 ) with b = Sj t ji. We can now conclude by writing down the 
full equations of motion for the coefficients for the two blob moments as: 



dAP[fci,fc 2 ;t] 
dt 



A(k 1 ,k 2 ) + B(k 1 ,k 2 ) + C(k 1 ,k 2 ), 



(A.35) 
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where 



(_l)(fci+te)^2(fci+fe 2 ) ~ ™ 
2 fci+fc2(fc 1 !)(fc 2 !) 



A{kiM)= v n '^ fl ',' ut „ E £ M^^i.^.tlAf^fmi.ma.^flfc^fca.^^a.mLmajA], (A.36) 



where r[fci, fc 2) 4, £ 2 , "H, m 2 ; A] is defined in equations (A. 19) and (A. 20). B(k\,k 2 ) represents the j 7^ j' 
terms in the first sum and can be written as 



B(k u k 2 ) = 



(_l)(fel+fe2)^2(fe 1 +fe 2 ) "\ °° 



2**+*»(jfc 1 !)(jfe a !) 



[mi,m 2 ,t]r B [h, k 2) £i,£ 2 ,mi,m 2 ; A, {sjj/ JA.37) 



j'^j £1/2 "11, "12 



where r B [Jfei , fc 2 ; 4 , 4 , mi , m 2 ; A, {s jd > }} = T 1 B [ki,k 2 ;£i,£ 2 ,mi,m 2 ;X, {s jtj i }} +rg [fc x , & 2 ; 4 , l 2 , m x , m 2 ; A, {s.y/ }] 
and 



min(£i ,ki — 1) min(£2-&2) 

E E 

i=Q j=0 

•Hf (mi + fci - i - 1 +4 - i,m 2 + k 2 -j +l 2 - j) 



, , i 2 i+1 fci! \ / 2 7,-j! 

>/\j) 2 \ A 2 ( l+1 )(4 — i — 1)! J VA 2(j) (fc2 - .?)! 



(A.38) 



r ii 

1 B 



and 



2^1! 



2J +1 fc>! 



mm(£i,fei) min(£2-,k2 — 1) 

§ § V*AjV V ^ VA 2 W(^'i -^7 VA 2(3+1) (fc2 -j-l)!, 

Z — u j — 

Hi{m 1 + k 1 -i + e 1 -i,m 2 + k 2 -j-l + £ 2 - j), (A.39) 



,B _ no, n a 2 J_ (~ & 2,^l) A _ 
' " & 1 ^ 2tT IbP 1 



The final term, C, in equation (A. 35 1 is simply the contribution from (A. 33 1 i.e. 

C = i J 1 (<)AP[fc 1 - l,k 2 ,t] +x J 2 (t)M j [k 1 ,k 2 - l,t}. 



(A.40) 



(A.41) 



The simplified equations of motion for the centers can also be written compactly as: 
dx^ 



m oo oo 



df 



= mJq.^ EE E M^"[4,4;i]M^[mi,m 2 ;t] 

L ' ' J j'=14A mi, 



j'=i Zi,t2 mi, ma 
\mi+m 2 "T/B / 



x( _ ir+m2r(mi +4 >m2 +4.) 



(A.42) 



Thus equations ( A. 35 )-( A.42 ) represent the reduced equations of motion for vorticity as given by the model. 
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